home *** CD-ROM | disk | FTP | other *** search
- // This may look like C code, but it is really -*- C++ -*-
- /*
- ************************************************************************
- *
- * Linear Algebra Package
- *
- * Basic Linear Algebra operations, level 1 & 2
- * concerning specifically vectors
- *
- * The present file is concerned with the operations which either
- * - specifically defined for vectors, such as norms
- * - some BLAS 1 & 2 operations that can be implemented more
- * efficiently than general operations on n*1 matrices
- *
- ************************************************************************
- */
-
- #include "LinAlg.h"
- #include <math.h>
- #include <builtin.h>
-
- #pragma implementation
-
- /*
- *------------------------------------------------------------------------
- * Specific vector constructors
- */
-
- #include <stdarg.h>
- // Make a vector and assign initial values
- // Argument list should contain DOUBLE values
- // to assign to vector elements. The list must
- // be terminated by the string "END"
- // Example: Vector foo(1,3,0.0,1.0,1.5,"END");
- Vector::Vector(const short lwb, const short upb, double iv1, ... )
- : Matrix(lwb,upb,1,1)
- {
- va_list args;
- va_start(args,iv1); // Init 'args' to the beginning of
- // the variable length list of args
- register int i;
- (*this)(lwb) = iv1;
- for(i=lwb+1; i<=upb; i++)
- (*this)(i) = (double)va_arg(args,double);
-
- assure( strcmp((char *)va_arg(args,char *),"END") == 0,
- "Vector: argument list must be terminated by \"END\" ");
- }
-
- /*
- *------------------------------------------------------------------------
- * Computing Vector norms
- */
-
- // Compute the 1-norm of the vector
- // SUM{ |v[i]| }
- double Vector::norm_1(void) const
- {
- is_valid();
-
- register double norm = 0;
- register REAL * vp;
-
- for(vp = elements; vp < elements + nelems; )
- norm += ::abs(*vp++);
-
- return norm;
- }
-
- // Compute the square of the 2-norm
- // SUM{ v[i]^2 }
- double Vector::norm_2_sqr(void) const
- {
- is_valid();
-
- register double norm = 0;
- register REAL * vp;
-
- for(vp = elements; vp < elements + nelems; )
- norm += ::sqr(*vp++);
-
- return norm;
- }
- // Compute the infinity-norm of the vector
- // MAX{ |v[i]| }
- double Vector::norm_inf(void) const
- {
- is_valid();
-
- register double norm = 0;
- register REAL * vp;
-
- for(vp = elements; vp < elements + nelems; )
- norm = ::max(norm,::abs(*vp++));
-
- return norm;
- }
-
- /*
- *------------------------------------------------------------------------
- * Multiplications specifically defined for vectors
- */
-
- // Compute the scalar product
- double operator * (const Vector& v1, const Vector& v2)
- {
- are_compatible(v1,v2);
- register REAL * v1p = v1.elements;
- register REAL * v2p = v2.elements;
- register double sum = 0;
-
- while( v1p < v1.elements + v1.nelems )
- sum += *v1p++ * *v2p++;
-
- return sum;
- }
-